//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Octant.cc
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Octant member definitions and template instantiations.
 * \note   Copyright (C) 2011 Jeremy Roberts.  Derived from Denovo.
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <cmath>

#include "angle/Octant.hh"
#include "tran/Definitions.hh"
#include "../utilities/Constants.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/DBC.hh"


namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Build an octant.
 *
 * A quadrature set is provided by the quadrature object.  A
 * valid Octant_Designator determines the octant.
 *
 * \param octant_in Octant_Designator that determines the octant
 * \param quadrature quadrature set
 *
 * \param do_adjoint if true, set the angles in each octant to the negative of
 * the octant direction in order to do adjoint sweeps
 */
template <class D>
Octant<D>::Octant(const Octant_Designator &octant,
                  const Quadrature<D>     &quadrature,
                  bool                     do_adjoint)
    : d_octant(octant)
    , d_dimension(D::dimension)
    , d_number_angles( quadrature.number_angles() / std::pow(2,d_dimension) )
    , d_angles(d_number_angles)
    , d_weights(d_number_angles, 0.0)
    , d_adjoint(do_adjoint)
{
    using util::soft_equiv;
    
    int num_octs = std::pow(2,d_dimension);
    Require (quadrature.number_angles() % num_octs == 0);
    for ( int i = 0; i < d_dimension; i++ )
    	Require (d_octant[i] == 1 || d_octant[i] == -1);

    // switch sign of the octant if in mode adjoint mode
    double adj = 1.0;
    if (d_adjoint) adj = -1.0;

    // load the quadrature angles
    load_quadrature_angles(quadrature,adj);

}

//! definition
template <class D>
void Octant<D>::load_quadrature_angles(const Quadrature<D> &quadrature, double adj)
{
	for ( int i = 0; i < d_number_angles; i++ )
	{
		// i indexes the Omega_i, composed of mu/eta/(xi) indexed by j
		for ( int j = 0; j < d_dimension; j++ )
			d_angles[i][j]  = adj * quadrature.angle(j)[i] * d_octant[j];
		d_weights[i] = quadrature.weight(i);
	}
}

//! specialization for 1D
template <>
void Octant<_1D>::load_quadrature_angles(const Quadrature<_1D> &quadrature, double adj)
{
    for ( int i = 0; i < d_number_angles; i++ )
    {
        d_angles[i]  = adj * quadrature.angle(0)[i] * d_octant[0];
        d_weights[i] = quadrature.weight(i);
    }
}

// instantiations
template class Octant<_1D>;
template class Octant<_2D>;
template class Octant<_3D>;

} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of Octant.cc
//---------------------------------------------------------------------------//
